## Section 3.2.4.A of the pre-analysis plan: "Photo-Specific Measures
## (conditional on Likert-type or pairwise data) A. Distance between
## Posterior Mean Ranks"
##
## Kevin Quinn
## University of Michigan
##
## 7/15/2019
##

#library(TopKLists)
library(coda)

set.seed(2102983)


########################################################################
## START a1-b1 comparisons
########################################################################

## load MCMC output
load("MM3.2.a1.M1.out.Rda")
a1.M1.out <- a1.M1.out[sample(1:nrow(a1.M1.out), nrow(a1.M1.out),
                              replace=FALSE), ]
load("MM3.2.a1.M2.out.Rda")
a1.M2.out <- a1.M2.out[sample(1:nrow(a1.M2.out), nrow(a1.M2.out),
                              replace=FALSE), ]
load("MM3.2.a1.pairwise.out.Rda")
a1.pairwise.out <- a1.pairwise.out[sample(1:nrow(a1.pairwise.out),
                                          nrow(a1.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b1.M1.out.Rda")
b1.M1.out <- b1.M1.out[sample(1:nrow(b1.M1.out), nrow(b1.M1.out),
                              replace=FALSE), ]

load("MM3.2.b1.M2.out.Rda")
b1.M2.out <- b1.M2.out[sample(1:nrow(b1.M2.out), nrow(b1.M2.out),
                              replace=FALSE), ]
load("MM3.2.b1.pairwise.out.Rda")
b1.pairwise.out <- b1.pairwise.out[sample(1:nrow(b1.pairwise.out),
                                          nrow(b1.pairwise.out),
                              replace=FALSE), ]

## make MCMC output comparable
a1.M1.out <- a1.M1.out[,1:48]
a1.M2.out <- a1.M2.out[,1:48]
b1.M1.out <- b1.M1.out[,1:48]
b1.M2.out <- b1.M2.out[,1:48]

colnames(a1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M1.out))
colnames(a1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M2.out))
colnames(b1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M1.out))
colnames(b1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M2.out))
colnames(a1.pairwise.out) <- gsub("theta\\.", "", colnames(a1.pairwise.out))
colnames(b1.pairwise.out) <- gsub("theta\\.", "", colnames(b1.pairwise.out))

for (j in 2:48){
    a1.M1.out[,j] <- a1.M1.out[,1] + a1.M1.out[,j]
    a1.M2.out[,j] <- a1.M2.out[,1] + a1.M2.out[,j]
    b1.M1.out[,j] <- b1.M1.out[,1] + b1.M1.out[,j]
    b1.M2.out[,j] <- b1.M2.out[,1] + b1.M2.out[,j]
}
colnames(a1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(a1.M2.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(b1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(b1.M2.out)[1] <- colnames(a1.pairwise.out)[1]


M <- nrow(a1.M1.out)
J <- ncol(a1.M1.out)

rank.a1.M1 <- matrix(NA, M, J)
rank.a1.M2 <- matrix(NA, M, J)
rank.b1.M1 <- matrix(NA, M, J)
rank.b1.M2 <- matrix(NA, M, J)
rank.a1.pair <- matrix(NA, M, J)
rank.b1.pair <- matrix(NA, M, J)

for (iter in 1:M){
    rank.a1.M1[iter,] <- rank(a1.M1.out[iter,])
    rank.a1.M2[iter,] <- rank(a1.M2.out[iter,])
    rank.b1.M1[iter,] <- rank(b1.M1.out[iter,])
    rank.b1.M2[iter,] <- rank(b1.M2.out[iter,])
    rank.a1.pair[iter,] <- rank(a1.pairwise.out[iter,])
    rank.b1.pair[iter,] <- rank(b1.pairwise.out[iter,])
}

diff.a1b1.M1 <- colMeans(rank.a1.M1) - colMeans(rank.b1.M1)
diff.a1b1.M2 <- colMeans(rank.a1.M2) - colMeans(rank.b1.M2)
diff.a1b1.pair <- colMeans(rank.a1.pair) - colMeans(rank.b1.pair)
ordvec <- order(diff.a1b1.pair)

pdf(file="MM3.2.4.A.a1b1.pdf", height=6, width=7)
plot(diff.a1b1.M1[ordvec], xlab="Photo Index",
     ylab="Difference in Expected Rank",
     main="Coder Race Correlates with Photos Observed", pch=0,
     col="orangered1", ylim=c(-17,17))
abline(h=0)
points(diff.a1b1.M1[ordvec], pch=0, col="orangered1")
points(diff.a1b1.M2[ordvec], pch=1, col="dodgerblue")
points(diff.a1b1.pair[ordvec], pch=3, col="purple4")
legend(x=5, y=17, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
       col=c("orangered1", "dodgerblue", "purple4"))
dev.off()
########################################################################
## END a1-b1 comparisons
########################################################################






########################################################################
## START a2-b2 comparisons
########################################################################

## load MCMC output
load("MM3.2.a2.M1.out.Rda")
a2.M1.out <- a2.M1.out[sample(1:nrow(a2.M1.out), nrow(a2.M1.out),
                              replace=FALSE), ]
load("MM3.2.a2.M2.out.Rda")
a2.M2.out <- a2.M2.out[sample(1:nrow(a2.M2.out), nrow(a2.M2.out),
                              replace=FALSE), ]
load("MM3.2.a2.pairwise.out.Rda")
a2.pairwise.out <- a2.pairwise.out[sample(1:nrow(a2.pairwise.out),
                                          nrow(a2.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b2.M1.out.Rda")
b2.M1.out <- b2.M1.out[sample(1:nrow(b2.M1.out), nrow(b2.M1.out),
                              replace=FALSE), ]

load("MM3.2.b2.M2.out.Rda")
b2.M2.out <- b2.M2.out[sample(1:nrow(b2.M2.out), nrow(b2.M2.out),
                              replace=FALSE), ]
load("MM3.2.b2.pairwise.out.Rda")
b2.pairwise.out <- b2.pairwise.out[sample(1:nrow(b2.pairwise.out),
                                          nrow(b2.pairwise.out),
                              replace=FALSE), ]

## make MCMC output comparable
a2.M1.out <- a2.M1.out[,1:48]
a2.M2.out <- a2.M2.out[,1:48]
b2.M1.out <- b2.M1.out[,1:48]
b2.M2.out <- b2.M2.out[,1:48]

colnames(a2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M1.out))
colnames(a2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M2.out))
colnames(b2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M1.out))
colnames(b2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M2.out))
colnames(a2.pairwise.out) <- gsub("theta\\.", "", colnames(a2.pairwise.out))
colnames(b2.pairwise.out) <- gsub("theta\\.", "", colnames(b2.pairwise.out))

for (j in 2:48){
    a2.M1.out[,j] <- a2.M1.out[,1] + a2.M1.out[,j]
    a2.M2.out[,j] <- a2.M2.out[,1] + a2.M2.out[,j]
    b2.M1.out[,j] <- b2.M1.out[,1] + b2.M1.out[,j]
    b2.M2.out[,j] <- b2.M2.out[,1] + b2.M2.out[,j]
}
colnames(a2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(a2.M2.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(b2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(b2.M2.out)[1] <- colnames(a2.pairwise.out)[1]


M <- nrow(a2.M1.out)
J <- ncol(a1.M1.out)

rank.a2.M1 <- matrix(NA, M, J)
rank.a2.M2 <- matrix(NA, M, J)
rank.b2.M1 <- matrix(NA, M, J)
rank.b2.M2 <- matrix(NA, M, J)
rank.a2.pair <- matrix(NA, M, J)
rank.b2.pair <- matrix(NA, M, J)

for (iter in 1:M){
    rank.a2.M1[iter,] <- rank(a2.M1.out[iter,])
    rank.a2.M2[iter,] <- rank(a2.M2.out[iter,])
    rank.b2.M1[iter,] <- rank(b2.M1.out[iter,])
    rank.b2.M2[iter,] <- rank(b2.M2.out[iter,])
    rank.a2.pair[iter,] <- rank(a2.pairwise.out[iter,])
    rank.b2.pair[iter,] <- rank(b2.pairwise.out[iter,])
}

diff.a2b2.M1 <- colMeans(rank.a2.M1) - colMeans(rank.b2.M1)
diff.a2b2.M2 <- colMeans(rank.a2.M2) - colMeans(rank.b2.M2)
diff.a2b2.pair <- colMeans(rank.a2.pair) - colMeans(rank.b2.pair)
ordvec <- order(diff.a2b2.pair)

pdf(file="MM3.2.4.A.a2b2.pdf", height=6, width=7)
plot(diff.a2b2.M1[ordvec], xlab="Photo Index",
     ylab="Difference in Expected Rank",
     main="Distractor-Photo Race Correlates with Photos Observed", pch=0,
     col="orangered1", ylim=c(-17,17))
abline(h=0)
points(diff.a2b2.M1[ordvec], pch=0, col="orangered1")
points(diff.a2b2.M2[ordvec], pch=1, col="dodgerblue")
points(diff.a2b2.pair[ordvec], pch=3, col="purple4")
legend(x=5, y=17, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
       col=c("orangered1", "dodgerblue", "purple4"))
dev.off()




########################################################################
## END a2-b2 comparisons
########################################################################







########################################################################
## START a3-b3 comparisons
########################################################################

## load MCMC output
load("MM3.2.a3.M1.out.Rda")
a3.M1.out <- a3.M1.out[sample(1:nrow(a3.M1.out), nrow(a3.M1.out),
                              replace=FALSE), ]
load("MM3.2.a3.M2.out.Rda")
a3.M2.out <- a3.M2.out[sample(1:nrow(a3.M2.out), nrow(a3.M2.out),
                              replace=FALSE), ]
load("MM3.2.a3.pairwise.out.Rda")
a3.pairwise.out <- a3.pairwise.out[sample(1:nrow(a3.pairwise.out),
                                          nrow(a3.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b3.M1.out.Rda")
b3.M1.out <- b3.M1.out[sample(1:nrow(b3.M1.out), nrow(b3.M1.out),
                              replace=FALSE), ]

load("MM3.2.b3.M2.out.Rda")
b3.M2.out <- b3.M2.out[sample(1:nrow(b3.M2.out), nrow(b3.M2.out),
                              replace=FALSE), ]
load("MM3.2.b3.pairwise.out.Rda")
b3.pairwise.out <- b3.pairwise.out[sample(1:nrow(b3.pairwise.out),
                                          nrow(b3.pairwise.out),
                              replace=FALSE), ]

## make MCMC output comparable
a3.M1.out <- a3.M1.out[,1:48]
a3.M2.out <- a3.M2.out[,1:48]
b3.M1.out <- b3.M1.out[,1:48]
b3.M2.out <- b3.M2.out[,1:48]

colnames(a3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M1.out))
colnames(a3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M2.out))
colnames(b3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M1.out))
colnames(b3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M2.out))
colnames(a3.pairwise.out) <- gsub("theta\\.", "", colnames(a3.pairwise.out))
colnames(b3.pairwise.out) <- gsub("theta\\.", "", colnames(b3.pairwise.out))

for (j in 2:48){
    a3.M1.out[,j] <- a3.M1.out[,1] + a3.M1.out[,j]
    a3.M2.out[,j] <- a3.M2.out[,1] + a3.M2.out[,j]
    b3.M1.out[,j] <- b3.M1.out[,1] + b3.M1.out[,j]
    b3.M2.out[,j] <- b3.M2.out[,1] + b3.M2.out[,j]
}
colnames(a3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(a3.M2.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(b3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(b3.M2.out)[1] <- colnames(a3.pairwise.out)[1]


M <- nrow(a3.M1.out)
J <- ncol(a1.M1.out)

rank.a3.M1 <- matrix(NA, M, J)
rank.a3.M2 <- matrix(NA, M, J)
rank.b3.M1 <- matrix(NA, M, J)
rank.b3.M2 <- matrix(NA, M, J)
rank.a3.pair <- matrix(NA, M, J)
rank.b3.pair <- matrix(NA, M, J)

for (iter in 1:M){
    rank.a3.M1[iter,] <- rank(a3.M1.out[iter,])
    rank.a3.M2[iter,] <- rank(a3.M2.out[iter,])
    rank.b3.M1[iter,] <- rank(b3.M1.out[iter,])
    rank.b3.M2[iter,] <- rank(b3.M2.out[iter,])
    rank.a3.pair[iter,] <- rank(a3.pairwise.out[iter,])
    rank.b3.pair[iter,] <- rank(b3.pairwise.out[iter,])
}

diff.a3b3.M1 <- colMeans(rank.a3.M1) - colMeans(rank.b3.M1)
diff.a3b3.M2 <- colMeans(rank.a3.M2) - colMeans(rank.b3.M2)
diff.a3b3.pair <- colMeans(rank.a3.pair) - colMeans(rank.b3.pair)
ordvec <- order(diff.a3b3.pair)

pdf(file="MM3.2.4.A.a3b3.pdf", height=6, width=7)
plot(diff.a3b3.M1[ordvec], xlab="Photo Index",
     ylab="Difference in Expected Rank",
     main="Black vs. White Coders", pch=0,
     col="orangered1", ylim=c(-17,17))
abline(h=0)
points(diff.a3b3.M1[ordvec], pch=0, col="orangered1")
points(diff.a3b3.M2[ordvec], pch=1, col="dodgerblue")
points(diff.a3b3.pair[ordvec], pch=3, col="purple4")
legend(x=5, y=17, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
       col=c("orangered1", "dodgerblue", "purple4"))
dev.off()




########################################################################
## END a3-b3 comparisons
########################################################################





########################################################################
## START a4-b4 comparisons
########################################################################

## load MCMC output
load("MM3.2.a4.M1.out.Rda")
a4.M1.out <- a4.M1.out[sample(1:nrow(a4.M1.out), nrow(a4.M1.out),
                              replace=FALSE), ]
load("MM3.2.a4.M2.out.Rda")
a4.M2.out <- a4.M2.out[sample(1:nrow(a4.M2.out), nrow(a4.M2.out),
                              replace=FALSE), ]
load("MM3.2.a4.pairwise.out.Rda")
a4.pairwise.out <- a4.pairwise.out[sample(1:nrow(a4.pairwise.out),
                                          nrow(a4.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b4.M1.out.Rda")
b4.M1.out <- b4.M1.out[sample(1:nrow(b4.M1.out), nrow(b4.M1.out),
                              replace=FALSE), ]

load("MM3.2.b4.M2.out.Rda")
b4.M2.out <- b4.M2.out[sample(1:nrow(b4.M2.out), nrow(b4.M2.out),
                              replace=FALSE), ]
load("MM3.2.b4.pairwise.out.Rda")
b4.pairwise.out <- b4.pairwise.out[sample(1:nrow(b4.pairwise.out),
                                          nrow(b4.pairwise.out),
                              replace=FALSE), ]

## make MCMC output comparable
a4.M1.out <- a4.M1.out[,1:48]
a4.M2.out <- a4.M2.out[,1:48]
b4.M1.out <- b4.M1.out[,1:48]
b4.M2.out <- b4.M2.out[,1:48]

colnames(a4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M1.out))
colnames(a4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M2.out))
colnames(b4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M1.out))
colnames(b4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M2.out))
colnames(a4.pairwise.out) <- gsub("theta\\.", "", colnames(a4.pairwise.out))
colnames(b4.pairwise.out) <- gsub("theta\\.", "", colnames(b4.pairwise.out))

for (j in 2:48){
    a4.M1.out[,j] <- a4.M1.out[,1] + a4.M1.out[,j]
    a4.M2.out[,j] <- a4.M2.out[,1] + a4.M2.out[,j]
    b4.M1.out[,j] <- b4.M1.out[,1] + b4.M1.out[,j]
    b4.M2.out[,j] <- b4.M2.out[,1] + b4.M2.out[,j]
}
colnames(a4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(a4.M2.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(b4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(b4.M2.out)[1] <- colnames(a4.pairwise.out)[1]


M <- nrow(a4.M1.out)
J <- ncol(a1.M1.out)

rank.a4.M1 <- matrix(NA, M, J)
rank.a4.M2 <- matrix(NA, M, J)
rank.b4.M1 <- matrix(NA, M, J)
rank.b4.M2 <- matrix(NA, M, J)
rank.a4.pair <- matrix(NA, M, J)
rank.b4.pair <- matrix(NA, M, J)

for (iter in 1:M){
    rank.a4.M1[iter,] <- rank(a4.M1.out[iter,])
    rank.a4.M2[iter,] <- rank(a4.M2.out[iter,])
    rank.b4.M1[iter,] <- rank(b4.M1.out[iter,])
    rank.b4.M2[iter,] <- rank(b4.M2.out[iter,])
    rank.a4.pair[iter,] <- rank(a4.pairwise.out[iter,])
    rank.b4.pair[iter,] <- rank(b4.pairwise.out[iter,])
}

diff.a4b4.M1 <- colMeans(rank.a4.M1) - colMeans(rank.b4.M1)
diff.a4b4.M2 <- colMeans(rank.a4.M2) - colMeans(rank.b4.M2)
diff.a4b4.pair <- colMeans(rank.a4.pair) - colMeans(rank.b4.pair)
ordvec <- order(diff.a4b4.pair)

pdf(file="MM3.2.4.A.a4b4.pdf", height=6, width=7)
plot(diff.a4b4.M1[ordvec], xlab="Photo Index",
     ylab="Difference in Expected Rank",
     main="Black vs. White Distractor Photos", pch=0,
     col="orangered1", ylim=c(-17,17))
abline(h=0)
points(diff.a4b4.M1[ordvec], pch=0, col="orangered1")
points(diff.a4b4.M2[ordvec], pch=1, col="dodgerblue")
points(diff.a4b4.pair[ordvec], pch=3, col="purple4")
legend(x=5, y=17, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
       col=c("orangered1", "dodgerblue", "purple4"))
dev.off()




########################################################################
## END a4-b4 comparisons
########################################################################

